#!/bin/env Rscript

library(HMMcopy)

argss <- commandArgs(trailingOnly = TRUE)

print(argss)


samplew <- argss[2]
sample_name <- argss[3]
gcw <- argss[4]
mapw <- argss[5]
chr <- argss[6]

chr=strsplit(chr, ',')[[1]]

# Correct readcount and estimate copy number
sample_uncorrected_reads <- wigsToRangedData(samplew, gcw, mapw)
# control_uncorrected_reads <- wigsToRangedData(controlw, gcw, mapw)

# The correctReadcount requires at least about 1000 bins to work properly
# process with chr separately?
# sample_corrected_copy <- correctReadcount(sample_uncorrected_reads)

source(paste(Sys.getenv('tools_path'), 'script/correct_rc_hmmcopy.r', sep='/'))
sample_corrected_copy <- correctrchmmcopy(sample_uncorrected_reads)
# control_corrected_copy <- correctReadcount(control_uncorrected_reads)

# Segmentation

# Normalizing sample by control
# sample_corrected_copy$copy <- sample_corrected_copy$copy - 1
# Segmenting

chr <- space(sample_corrected_copy)

autosomes <- (chr != "X" & chr != "Y" & chr != "23" & 
		chr != "24" & chr != "chrX" & chr != "chrY" & chr != 
		"M" & chr != "MT" & chr != "chrM")

# param <- data.frame(strength = 1e+07, e = 0.9999999, 
            # mu = quantile(correctOut$copy, na.rm = TRUE, prob = c(0.1, 
                # 0.25, 0.5, 0.75, 0.9, 0.99)), lambda = 20, nu = 2.1, 
            # kappa = c(0.05, 0.05, 0.7, 0.1, 0.05, 0.05) * 1000, 
            # m = 0, eta = c(5, 5, 50, 5, 5, 5) * 10000, gamma = 3, 
            # S = 0)
        # param$m <- param$mu
        # param$S <- ((sd(2^correctOut$copy[autosomes], na.rm = TRUE)/sqrt(nrow(param)))^2)
        # rownames(param) <- seq(1, 6)

mu = log(c(10^(-5), 1, 2, 3, 4, 5)/2, 2)
kappa = c(0.05, 0.1, 0.6, 0.15, 0.05, 0.05) * 1000
nu = 2.1
eta = c(5, 5, 50, 5, 5, 5) * 10000

param <- data.frame(strength = 1e+07, e = 0.9999999, mu, lambda = 20, nu, kappa, m = mu, eta, gamma = 3)

copy_range=quantile(sample_corrected_copy$copy[autosomes], c(0.05, 0.95), na.rm=T)

copy1=sample_corrected_copy$copy[sample_corrected_copy$copy > copy_range[1] & sample_corrected_copy$copy < copy_range[2]]

param$S <- ((sd(2^copy1, na.rm = TRUE)/sqrt(nrow(param)))^2)
        
sample_corrected_copy$copy[which(sample_corrected_copy$copy==-Inf)]=mu[1]
		
segmented_copy <- HMMsegment(sample_corrected_copy, param=param)
# retrieve converged parameters via EM
# param <- HMMsegment(sample_corrected_copy, getparam = TRUE)
# param$mu <- log(c(1, 1.4, 2, 2.7, 3, 4.5) / 2, 2)
# param$m <- param$mu
# perform segmentation via Viterbi
# segmented_copy <- HMMsegment(sample_corrected_copy, param)

# log2 states to normal states 
# set normal state to 2 
sample_corrected_copy$copy=2^(sample_corrected_copy$copy+1)
segmented_copy$segs$median=2^(segmented_copy$segs$median+1)


# Visualization

# chr <- c(1:22, "X", "Y")

source(paste(Sys.getenv('tools_path'), 'script/cnv_ideogram_p_hmmcopy.r', sep='/'))
cnvphmmcopy(sample_corrected_copy, segmented_copy, chr=chr, sample_name=sample_name)

q(save = "yes", status = 0, runLast = FALSE)
